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This paper presents a new approach to the estimation of the 
■ deformation of an isotropic Gaussian random field on R 2 based on 

dense observations of a single realization of the deformed random 
field. Under this framework we investigate the identification and es- 
timation of deformations. We then present a complete methodologi- 
f"H . cal package — from model assumptions to algorithmic recovery of the 

^0 ' deformation — for the class of nonstationary processes obtained by 

^ ,— < deforming isotropic Gaussian random fields. 

> 

1. Introduction. Random fields obtained by deforming the coordinates 
of an isotropic random field form a rich class of nonstationary processes. Such 
random fields have the form V(x) = Z(/ _1 (x)), where Z is a random field 
J> ' on M. d and / is a deformation, or a smooth bijection of M. d . In one dimension, 

for example, the deformation / is used in speech recognition to model local 
, compression or expansion of time in different utterances of a spoken word 

(see [17]). Working with deformations in more than one dimension, however, 
Tj- ■ has been a challenge. One dimensional deformations behave locally as a 

change of scale. In more than one dimension, however, deformations can 
rotate as well as scale local coordinates. These added dynamics, in addition 
to complicated requirements for invertibility, have restricted the pragmatic 
use of deformations when modeling processes in more than one dimension. 
^ ' In this paper we establish methodology for working with and estimating 
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deformations in M 2 from a single realization of a deformed isotropic Gaussian 
random field. 

The use of deformations to model nonstationary processes was first in- 
troduced to the spatial statistics literature by Sampson and Guttorp [18]. 
Their work, as well as that of subsequent authors (see, e.g., [7, 11, 16, 19]) 
deals mostly with data from sparse observation locations with independent 
replicates of the random field. Two recent papers by Clerc and Mallat [4, 5] 
consider a similar deformation model, but under a different observation sce- 
nario: a densely observed single realization. They use families of localized 
functions to estimate local properties of the deformation, and under the as- 
sumption of reflective shape recovery, these local properties are related to 
estimates of the shape of the reflective surface. Most of their results hold 
in one dimension, however, and it is not clear that in two dimensions their 
estimates work for arbitrarily smooth isotropic processes under general de- 
formations. 

Guyon and Perrin [10] tackle the problem of developing consistent esti- 
mates of deformations in two dimensions. They succeed in proving results 
within a class of deformations when observing random fields that are station- 
ary but not isotropic. In this paper we examine the consequence of adding 
the assumption of isotropy to the pre-deformed random fields. The isotropic 
assumption complicates the estimation of the original orientation of the lo- 
cal coordinates by introducing a rotational invariance to the random field. 
In Section 3 we notice the local behavior of C 1 diffeomorphisms can be ap- 
proximated by a composition of rotations, linear coordinate stretchings and 
translations. Estimating the initial local rotation — or the original orienta- 
tion of the local coordinates — using only observations in a local neighbor- 
hood becomes difficult. Our approach is to estimate two local parameters 
of the deformation, a dilatation and scale, which are invariant under initial 
local rotations. These parameters can be estimated from the local behavior 
of the deformed process and are still sufficient to uniquely characterize the 
deformation, up to global rigid motions. 

The structure of this paper is as follows. In Section 2 we present our 
modeling assumptions on the random fields to be deformed. We introduce a 
flexible semi-parametric class of random fields under which we can perform 
approximate likelihood estimation of the local parameters that character- 
ize the deformation. In Section 3 we make explicit our assumptions on the 
class of deformations and study in more detail the consequences of isotropy 
when estimating deformations. Sections 4 and 5 present the main contribu- 
tion of this paper. We outline our methodology and present an algorithm 
which allows fast construction of the deformation. Finally, we discuss some 
simulations and future work. 
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2. Deformation model. Our deformation model for nonstationarity is 
processes of the form Zo/" 1 , where Z is a constant mean isotropic Gaus- 
sian process on M 2 and / : M 2 — ► ]R 2 is an orientation preserving C 1 diffeo- 
morphism. We also add the assumption / has bounded local distortion, 
which implies our maps are quasiconformal (see Appendix 6 for the def- 
inition of quasiconformal and C 1 diffeomorphism) . By bounded local dis- 
tortion, we mean that the ratio of limsup x ^ xo |/(x) — /(xo)|/|x — xo| to 
liminfx-^xo l/( x ) ~~ /( x o)l/l x — x o| be uniformly bounded for all xo G M 2 . 

One of the advantages of this model is that estimates of f~ l can be used 
to transform the coordinates of Z o / , returning a nonstationary process 
to isotropy for which one can use existing statistical techniques. This can be 
done by mapping the observation locations Xj to / (xj), which transforms 
the graph of (x,;,Zo / _1 (xj)) to the graph of the original isotropic process 
Z at the locations / _1 (xj). 

We want the second order behavior of Z to be as general as possible while 
still allowing approximate likelihood techniques for estimation. We do this 
by introducing a regularization parameter a > which restricts the behavior 
of the autocovariance at the origin, and therefore controls the mean square 
smoothness of the isotropic processes Z. To simplify notation, let p a denote 
the nonnegative integer 



(1) Po 



[a/2], if a/2 ^Z, 
a/2-1, if a/2 6 Z. 



Let K be the autocovariance function for Z so that K (|t — s|) = cov{Z(t), Z(s)} 
for s, t G K 2 . We suppose there exists a constant c> and an a > so that 
JT(|t|) has 2p a continuous derivatives and 

(2) ^(l*D-E^^* 2fe ~ cG ^) «a|<|->0, 



where G a is defined by 
G a (t) -- 



(_l)i+K2j| t | Q; for a/2 ^Z, 

(_l)l+«/2| t |«i og |i| 5 for a/2 GZ. 

This class of autocovariances encompasses a broad range of processes in- 
cluding the Matern model (see [21]) and the so-called exponential family 
exp(— c 1 1 — s| 7 ) for 7 G (0,2). The parameter a controls the mean square 
differentiability of Z so that Z is n times mean square differentiable if and 
only if n < a/2. 

One advantage of this class of processes is that we can perform restricted 
maximum likelihood estimation under the distributional approximation of 
Z by an intrinsic random function of order [a/2\ with generalized co- 
variance function cG a (see [21]). In particular, suppose Z is a process 
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with autocovariance function K that satisfies (2) and suppose we observe 
Z := (Z(xi), . . . , Z(x n )). We can find a matrix L with n columns so that the 
covariance structure of LZ is easily approximated by a positive definite ma- 
trix that depends only on the function cG a . To find such a matrix L, we first 
need some notation: for x = (x,y) 6 R 2 and r = (?T,r2) G N 2 , let x r denote 
the monomial x ri y T ' 2 . Now let L be a matrix for which each row (u\, . . . ,u n ) 
satisfies J27=i u i' K i = f° r an r = ( r ii r 2) €E N such that r\ + r 2 < |_ a /2j- 
These rows can be easily computed by finding linearly independent vectors 
orthogonal to the space spanned by {(xj, . . . ,x£) : r± + V2 < [ a /2j }• Now let 
(tti, . . . , u n ) and (v±, . . . , v n ) be two rows of L and let k < [a/2\ , then 

n n 

^2^2u i v j \x i -x j \ 2k = 

i=lj=l 

so that 

cov ^UiZ(xi), 

U=i 

(3) 

n ii 

= c UiVjG a (\xi - Xj|) + Y UiVjR(\xi - Xj|), 

i,j = l i,j=l 

where R denotes the error term in the asymptotic approximation (2) so that 
= o(G a (\t\)) as |i| — > 0. Let (u{, . . . , u^) denote the £th row of L. Since 
the functions G a are conditionally positive definite of order La/2j, the ma- 
trix (cY^i t j=i ufu'jGadxi — Xj|)) P)9 is positive definite (see [21]). Therefore, 
by ignoring the last term in (3), we can approximate the covariance structure 
of LZ, under the semiparametric assumption (2) by the covariance matrix 
(cT,i,j=iUiU q j G a (\x i -x j \)) Ptq . 

We point out a useful fact that will be used in the following methodology: 
if one supposes the deformation / is smoother than the sample paths of the 
process Z, then Z and Z o f~ l share the same local mean squared smooth- 
ness. This is to our advantage in that a is invariant under deformations of 
the process, which should allow estimation of a directly from Z o f~ l with- 
out knowledge of /. For example, suppose the autocovariance of the process 
Z has the form (2) with a < 2 and suppose hGl 2 such that h / (0, 0). We 
show that there exists a constant < c\ < oo such that 

(4) E{Z(x + eh) -Z(x)} 2 - c x E{Z o /^(x + eh) - Z o f'H*)} 2 

as e — ► 0. In other words, the processes Z and Z o / _1 both have the same 
power decay of the variogram at the origin (see [6]), which can be estimated 
from one realization of Z o / . 

To see why (4) is true, first notice that 

(5) E{Z(x + eh) - Z(x)} 2 = 2K(0) - 2A'(|eh|) 
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as e — > (for the rest of this section we write f\ ~ fa to denote that fi/fa—*l 
as e — > with x and h fixed) . Similarly, 

E{Z o /-!( X + eh) - Z o /-'(x)} 2 ~ 2c|/- 1 (x + eh) - /^(x)] 01 

by direct application of (2) since (/^(x+eh)-/- 1 ^)! -^Oase^O. Now by 
the differentiability assumption on f~ l we also have |/ _1 (x+eh) — /~ 1 (x)| ~ 
C2€ for some < C2 < oo. This gives 

E{Z o /^(x + eh) - Z o / _1 (x)} 2 ~ 2 ccf e Q , 

which, in conjunction with (5), proves (4). 

The following methodology requires that the process Z be smoother than 
the deformation /. However, this assumption is made only to exclude difficul- 
ties when estimating a. In addition, we anticipate only using this method in 
practice for smooth deformations where the process Z is rarely very smooth. 
Therefore, we consider this assumption a minor technicality to the methods 
described below. 

Remark 1. In what follows we will be using approximate likelihood 
methods that will depend on the data only through sufficiently high order 
increments of the processes. Therefore, all our methods can be extended to 
intrinsic random functions (see [3]) such as fractional Brownian surfaces, for 
example. 

Remark 2. The dependence of K on the parameter c presents an iden- 
tifiability problem if one assumes the constant c is unknown. This is because 
c is confounded with changes of scale (which are included in our deformation 
class). For example, if the autocovariance for Z satisfies (2) with c = 1, then 
the autocovariance for the deformed processes Z(2x) satisfies (2) for c = 2 a . 
In what follows we choose to fix the constant c = 1 and assume it is known. 
This allows us to estimate the correct scale of the deformation. 

3. C 1 diffeomorphisms with bounded distortion. Since the original 
process, Z, is assumed to be isotropic, the most one can hope for is the 
identification of the deformation up to rigid motions of the plane. We will 
show in this section that this is indeed possible. We begin by studying the 
local behavior of the map / in terms of local affine transformations — whose 
existence is guaranteed by the assumption that / is a C 1 diffeomorphism. 
These local affine transformations can be decomposed using the singular 
value decomposition, which makes explicit how to obtain the local coordi- 
nates of the map by a composition of rotations and stretchings. Since the 
isotropic assumption on the original random field complicates the identifi- 
cation of the initial local rotation, we estimate the remaining parameters of 
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the local coordinates. These remaining parameters are characterized by an 
ellipse. The theory of quasiconformal maps, in part, studies maps through 
these ellipses, and provides us with the theoretical foundation we need. A 
brief overview of quasiconformal maps is included in Appendix 6. 

Since / is a diffeomorphism, its local linear behavior is characterized by 

(6) /(x + h) = /(x) + J / h + o(|h|), as|h|->0, 

where Jf is the Jacobian of / and, writing f(x,y) = (u(x,y),v(x,y)), is 
defined as 

m Ms ; 

Since / is orientation-preserving, det Jf > 0, and since / is C 1 , it has con- 
tinuous partial derivatives u x , u y , v x and v y . We can further decompose 
the map / by using the singular value decomposition to represent the linear 
map Jf as a sequence of rotations and stretchings. In particular, Jf = UAV t , 
where U, V are orthogonal matrices and A is a diagonal matrix with diag- 
onal elements Ai > A2 > (see [9]). Since det Jf > 0, we can take U, V to be 
rotation matrices, which gives 

cos 9 — sin#\ /Ai \ ( cosift — simft 



(8) j>=(; 



sin 9 cos 9 J \ A2 / \ sin tp cos ift 

where Ai > A2 > 0. This decomposition describes the local behavior of / by 
an initial rotation, a stretching and a final rotation. This sheds light on the 
information Z o f^ 1 provides about the map /. Since Z is rotationally invari- 
ant, the initial rotation by an angle of ift will be difficult, if not impossible, 
to estimate locally. However, our deformations map infinitesimal circles to 
infinitesimal ellipses and these ellipses are invariant under the initial local 
rotation of the map /, so there is some prospect of estimating them from 
one sample path of the process Z o The inclinations of these ellipses are 
given by 9 and their eccentricities by A1/A2 > 1. 

These locally defined ellipses are particularly useful in the study of de- 
formations. First, notice that the eccentricity provides a natural measure of 
local distortion that equals our original notion of distortion as the ratio of 
the upper and lower limits of |/(x) — /(x )|/|x — x | as x approaches x 
from different directions. This notion of distortion will become useful when 
studying the smoothing problem for local ellipse estimation. Second, there 
is existing literature that characterizes maps with a given ellipse field. For 
example, the inclinations and eccentricities of these ellipses are sufficient to 
recover f~ l up to postcomposition of conformal maps (see Appendix 6). 

We parameterize these ellipses by a pair {(//, eft) S C x M + : |/i| < 1}, where 
/i is referred to as the complex dilatation and (ft as the scale. Set \i to be the 
complex number 

(9) " = -^". 
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where K denotes the eccentricity of the ellipse and 6 denotes inclination. 
The scale parameter is defined by (f> = \/Ai^2- in the mathematical literature 
fi is called the complex dilatation of the map / . The complex dilatation of 
a map characterizes the infinitesimal ellipse that gets mapped to a infinites- 
imal circle. The mapping theorem from quasiconformal theory tells us that 
specifying (/i, <p) on a Jordan region is sufficient to characterize the map 
f^ 1 on Q up to rigid motions of Q (see Appendix 6). 

4. Methodology. This section outlines our procedure for the estimation 
of the deformation / _1 and the fractional index of the process a in the pres- 
ence of a densely observed sample path of a deformed Gaussian process. The 
estimation procedure will be done in three stages. First we estimate the frac- 
tional index of the process, taking advantage of the invariance of the fractal 
properties of the sample paths under sufficiently smooth deformations. Then, 
using the estimated fractional index d, we use a local likelihood approach to 
estimate the local complex dilatation and scale of the deformation, (/x,</>), 
independently over each neighborhood. Finally, we smooth and interpolate 
the estimates (/i, 4>) across neighborhoods. 

The two parameters (/i, 4>) are enough to characterize the deformation 
f^ 1 on the observation region up to a rotation and translation of the orig- 
inal coordinates. Therefore, on estimating (/i, (/>), one has estimated all the 
parameters in the model. To make this more useful, however, we need to 
recover / _1 from (^,<p). Sections 4.3.2 and 5 detail how one can efficiently 
reconstruct the deformation given (fJ,,(fi). 

For both the estimation of a and (/U, 4>), we use likelihood methods. Since 
we are supposing we have observed the process densely, any full likelihood 
method will be computationally impractical. Therefore, we partition the ob- 
servation locations into neighborhoods and assume independence of the pro- 
cess across partitions. To be clear, we are not changing our model assump- 
tions, just devising approximate likelihood techniques. Likelihood methods 
present two advantages for an initial study of estimation. First, there are 
no requirements on the configuration of observation locations. Second, the 
estimates are easily constructible and our experience has been that they are 
highly efficient. 

It is advantageous to switch to complex notation where a point (x,y) is 
represented by z = x + iy and consider our deformations as mapping regions 
of the complex plane. To set notation, we let z := (z\ , . . . , z n )' be the vector of 
observation locations in C, and Y := (Y±, . . . ,Y n )' denote the corresponding 
observations of the process Y := Z o f~ l . In what follows we dividing the 
observations into local neighborhoods. We denote by M n the partition of 
indices {1, . . . ,n} that corresponds to the local observation neighborhoods. 
For an index set I = {i\ , . . . , i m } G J\f n , let zj denote the vector (z^ , Zi m ) , 
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and similarly, define Yj. Last, for a function F(z,w), we let (F(z p , z q ))j 
denote the matrix with entries F(z p , z q ) for p,q £l. 

The description of the methodology is described as it applies to the sim- 
ulation pictured in Figure 1. It shows a simulation of a deformed process 
Y := Zo/ _1 , where the isotropic autocovariance function for Z has the form 
K (t) = 0.5151 - |t| a7 + 0(\t\ 2 ) as t -> 0; see [2] for details. The observations 
are on a grid of size 400 x 400 in [0, l] 2 and the deformation that transforms 
the coordinates of Y to isotropy is 

(10) f~\x + iy) := (1.2 - y) e"^ 1 "*)/ 2 + il.2. 

4.1. Estimating a. To estimate a, we first construct the neighborhood 
structure, Af n , of the observation locations. For this step alone, one would 
want the size of the neighborhoods to be as large as computationally feasible. 
Then a is obtained by maximizing an approximated log-likelihood under the 
supposition of independence across neighborhoods and that the process is 
isotropic with autocovariance satisfying (2). 

To approximate the log-likelihood for the observations Yj on a partic- 
ular neighborhood X = {i\, . . . ,i m } £ M n , we use the techniques discussed 
in Section 2 and approximate the log-likelihood for a linear transformation 
of the data LYj. We want to set the rows of L to be linearly independent 
vectors orthogonal to the space spanned by {{z\ , . . . , z\ m ) : r\ + T2 < [a/2j }, 
where (x + iy) r denotes the real number x ri y r2 . Unfortunately this involves 
knowing |_ck / 2 J . However, if we suppose we know a predetermined upper 
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Fig. 1. Deformed process. 
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bound on the magnitude of a, which we denote by a , it will be suf- 
ficient to find linearly independent vectors orthogonal to space spanned 
by {{z\ , . . . , z\ m ) : r\ + r2 < Loo/2j}- Now, using (3), the random vector 
Yj := LYj is multivariate Gaussian with covariance matrix 

EYjY^LiKQzp-z^V 

= L(G a (\z p - z g |)) x L' + L(R(\z p - z q \)) x L'. 

By ignoring the remainder term ~L(R(\z p — z q \))jli' and setting 
Sj jQ , := li(G a (\z p — z q \))j Jj', we get the following approximate log likelihood 
over neighborhood X: 

(11) e(a\Y x , z x ) = - i log |S x , a | - |Y x E^ Q Yi + C, 

where C is a constant additive term not depending on a. We suppose inde- 
pendence across blocks so the log-likelihood for the full observation vector Y 
is computed by summing the log-likelihoods (11) over neighborhoods 2 6 J\f . 
Now we estimate a by 

a:=argmax ^(a|Yj, zj). 

For the simulation in Figure 1, the neighborhood structure, M n , was con- 
structed by dividing the 400 x 400 grid into 40 x 40 blocks of size 10 x 10. 
The maximum approximate likelihood estimate is a = 0.696 (the true a is 
0.7). We do not address the accuracy of the estimates of a in this paper. 
However, in all our simulations we have seen very accurate estimation of a. 

4.2. Estimating (/jl,<P). In this section we suppose we have an estimate 
a so that by regarding this estimate as the truth, the only unknown param- 
eter is the quasiconformal map /, or equivalently, / . As was discussed in 
Section 3, we reparameterize / _1 by its complex dilatation and scale field 
(//, (j)). Fixing a, we use approximate maximum likelihood estimation on 
(//,</>) locally on each neighborhood to get a spatial map (£l,4>). 

Given a particular neighborhood I £ J\f n , we suppose the process behaves 
like a geometric anisotropic process Z(J^-ix) with smoothness parameter a. 
As was discussed in Section 3, only three parameters of Jt-i are estimable, 
and they can be reparameterized as a complex dilatation and scale (fii,(f>i). 
In complex notation the geometric anisotropic process can be written as 
Z(Ai(z + fJ^iz)), where \Ai\ 2 = cftf / (I — \fii\ 2 ) . Here i indexes the neighbor- 
hoods M n arid denotes these constants for each neighborhood. As in 
Section 4.1 consider the vector Yj := LYj, where the rows of L are linearly 
independent vectors orthogonal to {(z^ , . . . , z?) : r\ + < \aj2\ }. Letting 
9i denote (/Xj,^), we then estimate 6i by maximizing the log-likelihood 
^(0»|Yj,zx) = -log|Sx,eJ/2 - Y^E^.Yr/2 + C, where C is a constant 
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term not depending on 9{. Similar to what was done in Section 4.1, the 
covariance matrix Ex q. is computed by ignoring the remainder term so that 
£x,0, = \j{G & fi.{zp - z q ))xL', where G &i0i (z) := G & {\ A\ {z - Hiz)). Notice 
that the parameter range for fi is the complex unit disk {/x G C: \fi\ < 1} 
and 4> > 0. 

For the simulation in Figure 1, we again use the neighborhood structure 
of 40 x 40 blocks of size 10 x 10 to estimate The left-hand plot in 

Figure 2 graphs the local estimates fii as vectors. 

4.3. Smoothing and interpolating (fi,4>)- The parameters (fi,<f>) are lo- 
cally estimated values of the functions [i and (f>. Since the local likelihood 
estimation was done independently on each neighborhood, there was no 
smoothness constraint incorporated into the estimation of fi and (f>. To in- 
corporate smoothness conditions in these estimates, we choose to smooth 
(ft, <f>) after local likelihood estimation has been completed. As we will see, 
smoothing fx and (f> should, and will, be done in completely different ways. 

4.3.1. The complex dilatation field fi. The most naive approach to smooth- 
ing the dilatation field fi would be local averaging. However, it is worthwhile 
to investigate the smoothing problem under more generality. In particular, 
the dilatation field is a parameterization of the ellipse, which is the pre-image 
of a local infinitesimal circle under the deformation Just averaging the 
elements of a particular parameterization, however, does not provide a con- 
sistent notion of smoothing under reparameterization. In the following we 
develop a metric on the complex dilatation parameterization, based on a 
geometrical motivation, and use this metric to smooth ft. 
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Fig. 2. Le/t: 77ie estimated complex dilatation p,. Right: The result of smoothing the 
estimated dilatation fx. 
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The locally defined ellipses represent infinitesimal linear deformations, the 
magnitudes of which are naturally measured by their eccentricities, denoted 
by K. To develop a metric on ellipses, consider two matrix representations A 
and B of infinitesimal linear transformations. We claim that the magnitude 
of the distortion of the composed linear map AB^ 1 is a natural measure of 
distance between two ellipses. By (9), this is easily computed by 



-K-AB- 1 — ~, j 1 > 

1 - \Vab-A 



HA ~ HB 



1 - hah b 

Since K AB ~i is at least 1, we take logarithms to define the distance measure 
d on the complex dilatations: 

d(fi A , ^b) ■= \ \ogK AB -t = 1 log ( ] + Y AB - X \ , 
1 I \Y-\ii AB -i\J 

It interesting to note that d is indeed a metric and it equals the hyperbolic 
metric on the unit disk of C (see [12], e.g.). Much is known about this metric, 
in particular, the geodesic connecting two points is easily expressible as the 
arc on a circle that is orthogonal to the unit circle and joins the two points. 
One way to extend these ideas to define the convex combination, fj,*, of n 
complex dilatations with positive weights w\,. . . ,w n could be 

(12) n* =wgmm.{wid p (ii,Hi) H h w n dP(fi, fx n )} 

\n\<i 

for some power p > 0. Now by using weights that spatially vary, one can 
smooth the estimated dilatation field fi. In addition, convex combinations 
allow one to define bilinear or other interpolation approaches to produce 
maps of ft throughout the observation region Q. 

Unfortunately, we do not know of any closed form for (12), but a compu- 
tational approximation to the minimum is attainable. The right-hand plot 
of Figure 2 shows the result of the approximated minimization of (12) for 
p = 2 with uniform weights on a sliding window of size 4 x 4 in the interior, 
and the rectangular region that overlaps the sliding 4x4 window near the 
boundary. 

4.3.2. The scale field (p. Given a quasiconformal map /: — ► Q 1 with 
complex dilatation /i, all other quasiconformal maps with dilatation \x are 
obtained by postcomposing / with a conformal map on 0'. Therefore, once 
we smooth and interpolate jl, we need to use the scale information (j) to 
identify the deformation within the class of quasiconformal maps specified 
by /}. Section 5 outlines an algorithm for constructing a representative qua- 
siconformal map with prescribed dilatation. We use this map to transform 
the observation coordinates, then use the scale information to postcompose 
this map conformally to estimate the true deformation Figure 3 shows 
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an example of two maps with the same dilatation but different scale fields. 
The right-hand plot is the map / _1 given by equation (10) and the left-hand 
plot is a representative map, constructed from our algorithm, with the same 
dilatation /Uj-i. The two maps can be equated by postcomposing one with 
a conformal map. 

Consider deforming the observation space Q by a representative quasi- 
conformal map f:£l—>£l' with dilatation fj,* (see Section 5), where fj,* is 
the dilatation obtained by smoothing ft. Remember that fj,* is an estimate 
of fif-i so that / is an estimate — before postcomposition with a conformal 
map — of / . So by transforming the coordinates of the observed process 
Z o f^ 1 by /, the resulting process takes the form Zog -1 , where g = fof is 
a conformal map with local scale <fi g = \g'\ = <pf(4>f ° /)• Therefore, we want 
to find a conformal h that transforms g to the identity, that is, h = g~ l . 
Once we find such an h, the composition ho f will estimate f~ l up to rigid 
motions of C. In particular, we want | (h o g)'\ = \h' o g\ \g'\ = 1, so that 

(f>g°g 1 <t>f{<Pf°f) 

on Q'. Notice that our estimate 4> is a pointwise estimate of ^/-i = l/cfrfof^ 1 
at spatial locations z = (z±, . . . , z n ). Therefore, we want conformal h that 
satisfies 

(14) \h'( Wj )\ 



where (wi, ■ ■ ■ , w n ) := (f(zi),...,f(z n )) are points in 0'. Section 5 shows 
how one constructs / and computes <f)jof~ l efficiently. 
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Notice that h! is analytic on Q' and nonzero by invertibility. Therefore, 
log/i' is analytic on Q' , with real part log|/i'|. Since the real part of an an- 
alytic function is harmonic, by (14), we have that log<frj — logfij o f 
should approximate a sampled harmonic function at locations Wj . By rescal- 
ing, we suppose, without loss of generality, that £1 can be embedded into 
the unit disk D of the complex plane. Now, any harmonic function on the 
unit disk is harmonic when restricted to the region f2 C B. We then use the 
Fourier series to represent the boundary values for harmonic functions on 
B. Since the boundary values of a harmonic function uniquely determine its 
values on the interior, this gives a series representation of harmonic func- 
tions on £1. Unfortunately there are harmonic functions on Q that cannot 
be extended harmonically to the whole disk B. One way to avoid this diffi- 
culty is to conformally map f2 to the disk B instead of embedding it. This is 
guaranteed to work, but is more computationally expensive. For this reason, 
our examples used the embedding technique. 

We consider the class of harmonic functions on the unit disk represented 
by g{w) = Y^n=oA n w n for complex coefficients A n and some finite N . In po- 
lar coordinates (r,6), where w = re lS , this becomes g(r,6) = J2n=o A n r n e ind . 
We use this decomposition for the harmonic function log/i' on £1' . Since the 
real part of log b! is log | h'\ , we want to find a sequence A n such that 



N 



Re Kr"j n6 i ~ log 



n=0 <f>f°f ^Wj) 



for j = 1, . . . , n, where Wj = rje te j . Now by setting A n = a n + ib n , the problem 
becomes finding a n , b n so that 

N 1 



V r n Aa n cosnOj - b n sinnfy) w log J— — -. 

The advantage of this representation is the linearity in the parameters a n 
and b n . Using I 2 minimization to estimate a n + ib n , this corresponds to 
solving the linear problem min c ||Fc — ljU, where 

c := (ax, ...,a N ,bi,.. .,b N ), 
1 := flog 



<t>f°f 1 (wj)Jj 1 
F := [(r? cos n 9j)j m (-r^ sinn 9j)j n ]- 

The right-hand plot of Figure 4 shows the resulting smoothed (j) for the 
simulation in Figure 1. 

Once we have estimated the sequence A n , we can then approximate b! by 
h'(w) = exp{J2n=o A n w n }, which is then used to construct h. Finally, f~ l 
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Fig. 5. Left: Estimated deformation. Right: True deformation. 

is estimated by f~ l := ho f. Figure 5 plots f^ 1 (left) and the true defor- 
mation / . Now, using our estimated deformation, we can transform the 
coordinates of Y attempting to return the process to isotropy. The resulting 
process is shown in Figure 6. 

5. Generating a quasiconformal map from its complex dilatation. We 

first need to introduce some notation and present a few facts, most of which 
are explained in more detail in Appendix 6. The basic approach to our 
algorithm is to find a vector field flow representation of a quasiconformal 
map. By a vector field flow representation of a quasiconformal map /, we 
mean a class of maps indexed by a time variable, {ft}t€[o,i\i where /o is the 
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identity, f\ = f, and the functions ft depend smoothly on t so that 

(15) f t+£ o f t ~\z) = z + e {u t (z) + iv t (z)) + o(e) 

for some class of functions Ut, Vt- C — > R. By writing / = u + iv, Appendix 6 
shows that the complex dilatation of / can be expressed as the ratio djf /d z f, 



where d z f := ^(u x + v y ) + ^(v a 



and dzf :-- 



\ + %(v x + u y ). The 



final fact that we need for the presentation of the algorithm is an equation 
relating the complex dilatation of the composition g o / to the complex 
dilatations fj, g and fj,f. Using the chain rule to compute (h{g° f) an d d z (gof), 
we get 



l-lgof ■-- 

By rearranging terms, 
(16) 



&z(g o f) = Vf + (dzf/d z f)n g of 
d z {gof) 1 + (dj/d z f)n 9 o / ' 



fJ>g°f : 



dzf Vgof ~ Hf 



VgofVf 

In what follows we derive a set of differential equations for the vector fields 
{(iH, Vt)}te[o,i] which will be numerically solved to reconstruct /. A contin- 
uous extension of a construction in [1], page 99, demonstrates that one can, 
indeed, embed a quasiconformal map / into a vector field flow by specifying 
a class of dilatations {nt}t&[o,i] with boundary conditions fio = and [i\ = [x 
that depends smoothly on t. Let {ft}te[0,i] De a vector field flow such that 




9 2 0.4 0.0 OB 



Fig. 6. Estimated original isotropic process. 
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ft has dilatation fi t with condition fiQ = and \x\ = fi. Let {(ut, ^t)}je[o,i] 
denote the vector field associated with the flow ft so that equation (15) is 
satisfied. By setting g = ft+ £ o ff 1 and f = ft in the composition formula 
(16), one can easily compute the dilatation of ft+ £ ° ft' 1 as 

n -s ( Ht+e ~ IH d z ft\ ! / dtm d z f t \ , * 

Vi- m+entdzft/ \i-N 9*/*/ 

We can also compute the complex dilatation of z + e(u t + w 4 ) + o(e), by 
taking derivatives to give 

mo\ e{d x u t - dyv t ) + ie(d y u t + d x v t ) + o(e) . 
2 + e(9 x tit + ^wj) + ie{-d y ut + c^v*) + o(e) 

Then equating (17) with (18) and letting e — > gives 

dm_dzft_\ ^ t _x 
l-\Vt\ 2 dJt 
In particular, the vector field (ut,vt) satisfies 

(19) d x u t - d y v t = 2Re(a t ), 

(20) dyu t + d x v t = 2 Im(crt), 
where 

/ dm dj t \ j 

Now, given {fJ>t}te[o,i]-i ^ one can solve (19) and (20) for (ut,vt), then one 
can construct a time varying vector field flow realization of the map /. The 
usefulness of this representation is in its recursive nature. First notice that 
the initial map /o is the identity. Now supposing one has constructed ft 
up to some fixed time t < 1, one can compute tr f easily by deforming dtUti 
fit and d z ft/d z ft along with f t . Then by solving (19) and (20), one can 
approximate f t+£ by f t+£ = ft + e{u t ° ft + m o f t ) + o(e). 

To uncouple equations (19) and (20), represent (ut,Vt) by two functions 
<3? t and using u t = d y $ t + d x ^ t and v t = d x $ t - d y ^> t - Notice that ® t 
is the potential function and is the stream function for the vector field 
(vt,Ut). Equations (19) and (20) then become the Poisson equations 

(21) A*i = 2Re(<7 t ), 

(22) A* t =2Im(ff t ). 

To get a unique solution for these equations, one must specify boundary 
conditions on $t and *$>f These boundary conditions correspond to different 
vector field embeddings under the dilatation constraint jit- 
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The examples in the previous sections set fit '■= tfx for t G [0, 1] and the 
boundary conditions $t = *$>t = on dQti where Qt = /t(^)- In our imple- 
mentation of this algorithm, at each time step t, we interpolated at to a grid 
enclosing f^, then used a fast Poisson equation solver for grid data with 
Matlab. 

6. Discussion. Our main motivation for developing this methodology is 
to demonstrate that estimating deformations of arbitrarily smooth isotropic 
random fields is not only theoretically possible, but can be accomplished 
in practice. We make no claim to optimality of our methods, but we do 
hope that this will serve as a test case in the pursuit of efficient and robust 
methodology for these models in spatial statistics. The use of local likeli- 
hood techniques gives quick and easy estimates of (//, eft) that seem to work 
quite well in simulation. Unfortunately they are not very amenable to theo- 
retical study. We do think, however, that quasiconformal theory presents a 
promising avenue for the theoretical study and quantification of variability 
for estimates of deformations. In particular, we believe that quasiconformal 
theory will be useful for constructing consistent estimates of C 1 diffeomor- 
phisms and can also be used for developing flexible parametric models of C 1 
diffeomor phisms . 

We conclude this section with some simulations. First, we demonstrate 
our methods on differentiable processes. Figure 7 shows the results of apply- 
ing our methodology to a simulation of a deformed differentiable isotropic 
process. The true deformation is generated by a vector field flow and the 
original isotropic process has autocovariance function K of the form K(t) = 
0.0231 - (0.4034)t 2 + \t\ 3 + o{\t\ 3 ) as \t\ -> and is simulated on a 400 x 400 
square grid on [0, l] 2 ; see [2] for details. Notice that even though the process 
is smooth and, thus, it is more difficult to see the deformation, there is suf- 
ficient information embedded in higher order increments for very accurate 
estimation of the deformation. 

Another aspect of these models that warrants further study is the con- 
struction of distances on deformations that would allow comparison of method- 
ologies and a study of the variability of the deformation estimates. Since the 
deformations are only estimable up to rigid motions of the plane, we want 
these measures to be invariant under such postcomposition. Letting O de- 
note the observation region in C, the following are two candidates that have 
this invariance: 






Fig. 7. Deformation estimation for differ ■entiable processes. Upper left: The deformed 
process Z o f~ . Upper right: The estimate of the original isotropic process Z using the 
estimated deformation Bottom left: The true deformation f~ . Bottom right: The 

estimated deformation 

The first distance is based on how well the interpoint distances of / match 
/. The second distance is a function of the dilatations of / and /, instead 
of the pointwise values. 

As a first step toward quantifying the variability of the deformation esti- 
mates, we independently simulated the same deformed process as in Section 
4 three times. We then added white noise to the last two simulations with 
a standard deviation of 10% and 25% of the process standard deviation, 
respectively. The simulations and the corresponding deformation estimates 
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Table 1 

Two different distance measures on deformations (columns). The rows 
correspond to the simulations and estimates shown in Figure 8 





diC/- 1 ,/- 1 ) 




No noise 


0.0563 


0.0675 


10% noise 


0.0655 


0.0842 


25% noise 


0.1257 


0.1354 



are shown in Figure 8. First notice that the estimated deformation in Figure 
5 and the estimated deformation shown top right in Figure 8 are from the 
same estimation procedure applied to two independent samples. This gives a 
sense for the variability of the methodology applied to this deformed process. 
Notice the large scale structure of the deformation looks to be qualitatively 
stable. Also, notice that one could improve the deformation estimates from 
the last two samples by modeling the white noise in the likelihood estimation 
steps. We find it interesting that even without modeling the white noise the 
estimates are reasonable. 

For each estimate, we also computed its distance from the truth for each 
of our proposed distance measures (the integrals being approximated by 
Riemann sums). The results are displayed in Table 1. Notice that adding 
white noise degrades the deformation estimates. This can be seen both qual- 
itatively and from the distance measurements. We think, however, that the 
noisy estimates are reasonable and show a certain robustness and stability. 

APPENDIX: 

QUASICONFORMAL MAPS AND COMPLEX DILATATIONS 

The goal of this section is to give a short introduction to the theory of qua- 
siconformal maps and to state, without proof, the Mapping theorem (Ahlfors 
[1], Chapter 5) concerning the uniqueness of a map with a given dilatation. 
A complete study of quasiconformal maps can be found in Ahlfors [1], Lehto 
and Virtanen [15], Krushkal' [13] or Lawrynowicz [14]. The Mapping theo- 
rem will show that the complex dilatation, fj, of an orientation preserving 
C 1 diffeomorphism, /, with bounded distortion uniquely determines / up 
to postcomposition with a conformal map. The Mapping theorem is more 
general, however, and will not, for example, require the maps to be differen- 
tiable everywhere as in the case of C 1 diffeomorphisms. We first define C 1 
diffeomorphisms and take advantage of the extra smoothness assumptions 
to redefine /i. We then derive some properties of the complex dilatation 
and show how /z characterizes the infinitesimal ellipse that gets mapped to 
an infinitesimal circle under /. Finally, we will relax the C 1 diffeomorphic 
assumptions to state the Mapping theorem. 
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Fig. 8. Independent simulations of the deformed process from Section 4 with white noise 
(left column) and the corresponding estimated deformation (right column). 
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We will be considering maps / that are denned, not only on C, but also 
on domains C C. We choose to consider only well-behaved domains, which 
will allow succinct exposition. In particular, by a region, we mean a subset 
OcC that is either the whole plane C or is a simply connected open subset 
of C whose boundary, d£l, is a Jordan closed curve (also known as a Jordon 
region, see [20], e.g.). Define a homeomorphism / from a region f2 to £1' to 
be a continuous, one-to-one, onto map from Q to £1' such that both / and 
/ _1 are continuous. All the maps we wish to consider are necessarily home- 
omorphic; however, homeomorphisms are not smooth enough to allow us to 
talk about directional derivatives and Jacobians, which provide the easiest 
way to understand a complex dilatation. For this reason, we first restrict our 
homeomorphisms to be sufficiently smooth so that derivatives are defined 
everywhere and are continuous, that is, C 1 diffeomorphic. A homeomor- 
phism /(x) = (u(x), i>(x)) is an orientation preserving C 1 diffeomorphism if 
u and v are continuously differentiable and for any x£fl, 

(23) /(x + h) = /(x) + J f h + o(\h\) as|h|->0, 

where Jf is the Jacobian of the map / and det Jf > 0. Now since all the 
derivatives exist, we can define d z f and d%f as in Section 5. Notice that 
d z f 7^ since det J/ > and det Jf = u x v y — u y v x = \d z f\ 2 — \chf\ 2 - This 
allows us to define the complex dilatation // = ///:= c\f /d z f. We will see 
that n(z) agrees with the complex dilatation of the map defined in Section 3 
as parameterizing the eccentricity and inclination of the infinitesimal ellipse 
centered at z that gets mapped to an infinitesimal circle under /. 

To see how \i characterizes the local behavior of the diffeomorphism /, 
we switch to complex notation, z = x + iy, so we can write the behavior of 
/ as 

f( z ) = f( z o) + fx(zo)(x - X ) + fy{zo){y ~ 2/0 ) + o{\z - Z \). 

Rearranging terms, we get f(z) = f(zo) + d z f(z )(z-z ) + djf(z )(z - z ) + 
o{\z — zq\) and then factoring the nonzero term d z f, we can decompose the 
local behavior of / into an initial stretch map, a final rotation and a uniform 
stretching 

(24) f(z + h) = f(z ) + d z f(z )(h + fih) + o{\h\), 

where \i = d%f /d z f is the complex dilatation. Now define the directional 
derivative of / at an angle (3 £ [0, 2tt) by 

w = lim /(£±ff!!W(£>, 



the limit existing by the differentiability of /. Notice that when \x is zero 
at a point zq, (24) tells us that the directional derivative of / at zq does 
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not depend on direction, that is, that / is conformal at zq. In particular, 
dzf = reduces to the Cauchy-Riemann equations. Now since dpf(z) = 
e-^yiz + ee^)^, we get dpf = e"^/* cos 0+/,, sin/3) = 8J + e~ 2i ' 3 d z f . 
Therefore, the magnitudes of maximal and minimal stretching are attained 
at max^ \dpf\ = \d z f\ + \c\f\ and min^ \dpf\ = \d z f\ - \dgf\, where the last 
equality holds since det J/ = \d z f\ 2 — \c\f\ 2 > implies that \d z f\ > \&zf\. 
The singular value decomposition of Jj yields the representation (8) with 
Ai > A2 > 0, which equal the major and minor magnitudes of local stretching. 
Therefore, 

Ai _ max/? I dp f I _ 1 + 
A 2 min^ \dpf\ 1 — |a*| ' 

To see how \i is related to the inclination of the ellipse, notice that dpf = 
d z f{l + fie- 2i P) so that max^ \dpf\ is attained when 2/3 = arg(/x). Therefore, 
the inclination of the ellipse that gets mapped to the circle is arg(— /x)/2, 
which agrees with (9). 

In Section 3 we defined a measure of distortion at a point zq, induced by 
/, to be the ratio of limsup z ^ zo \ f(z) - f(z )\/\z - z \ to liminf 2 ^ \f(z) - 
f(zo)\/\z — zq\. Since this ratio is equal to max^ \dpf\/ minp \dpf\, our orig- 
inal notion of bounded distortion at a point zq is equivalent to |/i(^o)| < 1- 
Therefore, the assumption that / has uniformly bounded distortion amounts 
to the supposition that ||/^/||oo < 1- 

It turns out that our smoothness assumptions on / are unnatural and 
the presentation of the theory is best done in full generality. Instead of 
forcing the derivatives u x , u y , v x and v y to exist everywhere, we only suppose 
existence in the distributional sense and that they be locally integrable (see 
[8] ) . Here is the full definition of a quasiconformal map taken from [8] : 

Definition A.l. An orientation preserving homeomorphism f from a 
region Q to a region f(£l) is quasiconformal if there exists k < 1 such that 
f has locally integrable, distributional derivatives d z f and t\f on f2, and 
\c\f \ <k\d z f \ almost everywhere. 

In particular, any orientation preserving C 1 diffeomorphism with uni- 
formly bounded distortion is by definition considered a quasiconformal map. 
Now we have the following existence and uniqueness theorem, taken from 
[15], page 194. 

Theorem A.l. Let Q and Q' be conformally equivalent regions and [i a 
measurable function in O with ||//(^)||oo < 1- Then there exists a quasicon- 
formal mapping _/* : — ^ 0' whose complex dilatation coincides with fx almost 
everywhere. This mapping is uniquely determined up to a conformal mapping 
of O' onto itself. 
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A useful restatement of the above theorem is that all quasiconformal maps 
on Q can be represented by finding a map /, with dilatation /x, then post- 
composing it with conformal maps. 

We conclude by mentioning that for a general dilatation \x the associated 
quasiconformal maps will always be a homeomorphism, but not generally a 
C 1 diffeomorphism. There are sufficient conditions on the complex dilata- 
tion // to guarantee C 1 diffeomorphic solutions; see [15], Theorem 7.2, for 
example. 
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